Datasets
Models
Downloads: 1
884 lines (883 with data), 63.6 kB
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ResNet50 model \n",
"\n",
"**Due to GPU quota is only 30 hours/per week on Kaggle, each training need 15+ hours, so the notebook cann't commiting(otherwise will exceeding the quota), only download the csv files to submit**\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
"_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Using TensorFlow backend.\n"
]
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import pydicom\n",
"import os\n",
"import matplotlib.pyplot as plt\n",
"import collections\n",
"from tqdm import tqdm_notebook as tqdm\n",
"from datetime import datetime\n",
"\n",
"from math import ceil, floor, log\n",
"import cv2\n",
"\n",
"import tensorflow as tf\n",
"import keras\n",
"\n",
"import sys\n",
"\n",
"from keras_applications.resnet import ResNet50\n",
"\n",
"from sklearn.model_selection import ShuffleSplit\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model Parameters Setup"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"# Paths \n",
"Path = '../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/'\n",
"train_img_path = Path + 'stage_2_train/'\n",
"test_img_path = Path + 'stage_2_test/'\n",
"sample_csv = Path + \"stage_2_sample_submission.csv\"\n",
"train_csv = Path + 'stage_2_train.csv'\n"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"def read_testset(filename):\n",
" df = pd.read_csv(filename)\n",
" df[\"Image\"] = df[\"ID\"].str.slice(stop=12)\n",
" df[\"Diagnosis\"] = df[\"ID\"].str.slice(start=13)\n",
" \n",
" df = df.loc[:, [\"Label\", \"Diagnosis\", \"Image\"]]\n",
" df = df.set_index(['Image', 'Diagnosis']).unstack(level=-1)\n",
" \n",
" return df\n",
"\n",
"\n",
"\n",
"def read_trainset(filename = Path + \"stage_2_train.csv\"):\n",
" df = pd.read_csv(filename)\n",
" df[\"Image\"] = df[\"ID\"].str.slice(stop=12)\n",
" df[\"Diagnosis\"] = df[\"ID\"].str.slice(start=13)\n",
" duplicates_to_remove = [56346, 56347, 56348, 56349,\n",
" 56350, 56351, 1171830, 1171831,\n",
" 1171832, 1171833, 1171834, 1171835,\n",
" 3705312, 3705313, 3705314, 3705315,\n",
" 3705316, 3705317, 3842478, 3842479,\n",
" 3842480, 3842481, 3842482, 3842483 ]\n",
" df = df.drop(index = duplicates_to_remove)\n",
" df = df.reset_index(drop = True) \n",
" df = df.loc[:, [\"Label\", \"Diagnosis\", \"Image\"]]\n",
" df = df.set_index(['Image', 'Diagnosis']).unstack(level=-1)\n",
" return df\n",
"\n",
"\n",
" \n",
"test_df = read_testset(sample_csv)\n",
"train_df = read_trainset(train_csv)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead tr th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe thead tr:last-of-type th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr>\n",
" <th></th>\n",
" <th colspan=\"6\" halign=\"left\">Label</th>\n",
" </tr>\n",
" <tr>\n",
" <th>Diagnosis</th>\n",
" <th>any</th>\n",
" <th>epidural</th>\n",
" <th>intraparenchymal</th>\n",
" <th>intraventricular</th>\n",
" <th>subarachnoid</th>\n",
" <th>subdural</th>\n",
" </tr>\n",
" <tr>\n",
" <th>Image</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <td>ID_000000e27</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" </tr>\n",
" <tr>\n",
" <td>ID_000009146</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" </tr>\n",
" <tr>\n",
" <td>ID_00007b8cb</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" </tr>\n",
" <tr>\n",
" <td>ID_000134952</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" </tr>\n",
" <tr>\n",
" <td>ID_000176f2a</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" <td>0.5</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Label \\\n",
"Diagnosis any epidural intraparenchymal intraventricular subarachnoid \n",
"Image \n",
"ID_000000e27 0.5 0.5 0.5 0.5 0.5 \n",
"ID_000009146 0.5 0.5 0.5 0.5 0.5 \n",
"ID_00007b8cb 0.5 0.5 0.5 0.5 0.5 \n",
"ID_000134952 0.5 0.5 0.5 0.5 0.5 \n",
"ID_000176f2a 0.5 0.5 0.5 0.5 0.5 \n",
"\n",
" \n",
"Diagnosis subdural \n",
"Image \n",
"ID_000000e27 0.5 \n",
"ID_000009146 0.5 \n",
"ID_00007b8cb 0.5 \n",
"ID_000134952 0.5 \n",
"ID_000176f2a 0.5 "
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"test_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead tr th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe thead tr:last-of-type th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr>\n",
" <th></th>\n",
" <th colspan=\"6\" halign=\"left\">Label</th>\n",
" </tr>\n",
" <tr>\n",
" <th>Diagnosis</th>\n",
" <th>any</th>\n",
" <th>epidural</th>\n",
" <th>intraparenchymal</th>\n",
" <th>intraventricular</th>\n",
" <th>subarachnoid</th>\n",
" <th>subdural</th>\n",
" </tr>\n",
" <tr>\n",
" <th>Image</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <td>ID_000012eaf</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>ID_000039fa0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>ID_00005679d</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>ID_00008ce3c</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>ID_0000950d7</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Label \\\n",
"Diagnosis any epidural intraparenchymal intraventricular subarachnoid \n",
"Image \n",
"ID_000012eaf 0 0 0 0 0 \n",
"ID_000039fa0 0 0 0 0 0 \n",
"ID_00005679d 0 0 0 0 0 \n",
"ID_00008ce3c 0 0 0 0 0 \n",
"ID_0000950d7 0 0 0 0 0 \n",
"\n",
" \n",
"Diagnosis subdural \n",
"Image \n",
"ID_000012eaf 0 \n",
"ID_000039fa0 0 \n",
"ID_00005679d 0 \n",
"ID_00008ce3c 0 \n",
"ID_0000950d7 0 "
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"train_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Data EDA and Cleaning\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def correct_dcm(dcm):\n",
" x = dcm.pixel_array + 1000\n",
" px_mode = 4096\n",
" x[x>=px_mode] = x[x>=px_mode] - px_mode\n",
" dcm.PixelData = x.tobytes()\n",
" dcm.RescaleIntercept = -1000\n",
"\n",
"def window_image(dcm, window_center, window_width):\n",
" \n",
" if (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0) and (int(dcm.RescaleIntercept) > -1000):\n",
" correct_dcm(dcm)\n",
" \n",
" img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n",
" img_min = window_center - window_width // 2\n",
" img_max = window_center + window_width // 2\n",
" img = np.clip(img, img_min, img_max)\n",
"\n",
" return img\n",
"\n",
"def bsb_window(dcm):\n",
" brain_img = window_image(dcm, 40, 80)\n",
" subdural_img = window_image(dcm, 80, 200)\n",
" soft_img = window_image(dcm, 40, 380)\n",
" \n",
" brain_img = (brain_img - 0) / 80\n",
" subdural_img = (subdural_img - (-20)) / 200\n",
" soft_img = (soft_img - (-150)) / 380\n",
" bsb_img = np.array([brain_img, subdural_img, soft_img]).transpose(1,2,0)\n",
"\n",
" return bsb_img\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def window_with_correction(dcm, window_center, window_width):\n",
" if (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0) and (int(dcm.RescaleIntercept) > -100):\n",
" correct_dcm(dcm)\n",
" img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n",
" img_min = window_center - window_width // 2\n",
" img_max = window_center + window_width // 2\n",
" img = np.clip(img, img_min, img_max)\n",
" return img\n",
"\n",
"def window_without_correction(dcm, window_center, window_width):\n",
" img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n",
" img_min = window_center - window_width // 2\n",
" img_max = window_center + window_width // 2\n",
" img = np.clip(img, img_min, img_max)\n",
" return img\n",
"\n",
"def window_testing(img, window):\n",
" brain_img = window(img, 40, 80)\n",
" subdural_img = window(img, 80, 200)\n",
" soft_img = window(img, 40, 380)\n",
" \n",
" brain_img = (brain_img - 0) / 80\n",
" subdural_img = (subdural_img - (-20)) / 200\n",
" soft_img = (soft_img - (-150)) / 380\n",
" bsb_img = np.array([brain_img, subdural_img, soft_img]).transpose(1,2,0)\n",
"\n",
" return bsb_img\n"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADHCAYAAAAXg5iPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXl8VNXd/99n7sydPZlkspNAIBAEiaBRFAVFqSg+UK22WG3dqo+te5WqrbUtfWqtS221dW+1uLRYq9bt50ILVUEQMAqCQAKBQEL2SSaZfb2/P85NSCACUoIs9/16zSszd+49d5nPOTnLdxGapmFgYGBgcPhi+qovwMDAwMBgcDEaegMDA4PDHKOhNzAwMDjMMRp6AwMDg8Mco6E3MDAwOMwxGnoDAwODwxyjoT+IEEI8LoT42f7edw/llAohNCGE+b8ty8DgUEEIUSeE+NpXfR0HCqNyH0RomvaDwdjXwOBwQggxFXhe07Tir/paDhWMHv1BghBC+aqvwcBgMBhotGiMIA8sRkM/yAghxggh3hNC+IUQnwshvq5vnyeEeEwI8ZYQIgScrm+7q8+xtwkhmoQQjUKIq/QplpF9jr9Lfz9VCNEghJgjhGjVj7miTzn/I4T4VAjRLYSoF0LMPbBPweBQRghRIoR4RQjRJoTwCSEeFkKYhBB3CiG26pp7VgiRqe/fMx14pRBiG7BooG36vicJIZbq9WO13lvvOW+2EOIvuv47hRCvCiGcwNtAkRAiqL+K9Ov5sRCiVr/GF4UQ2X3KukS/Vp8Q4qcH9AEeBBgN/SAihLAAbwALgDzgBuCvQojR+i4XA78G3MCSnY49G7gF+BowEjhtD6crADKBIcCVwCNCiCz9uxBwKeAB/ge4Rghx3n91cwZHBPpI801gK1CK1NcLwOX663RgBOACHt7p8NOAMcBZA20TQgwB/h9wF5AN/Ah4WQiRq+/7HOAAjkbWn99rmhYCZgCNmqa59FcjcCNwnl5+EdAJPKLfw1jgMeAS/TsvcGRN+2iaZrwG6QVMAZoBU59t84G5wDzg2Z32nwfcpb9/GvhNn+9GAhowcoB9pwIRwNxn/1bgpC+4rgeRlQZk5dX6Hmu8jFfPC5gEtO2sD2AhcG2fz6OBBHLdr0dTI/p8P9C224Hndir3XeAyoBBIA1kDXNNUoGGnbeuBaX0+F/a5np8DL/T5zgnEga991c/3QL2MebLBpQio1zQt3WfbVmSvCKB+D8d+3Ofz7vYF8GmaluzzOYzsZSGEOBG4BxgHqIAV+Mcer97AAEqArTtpC6Q+t/b5vBXZqOb32TaQZvtuGwZ8Swgxq882C/Af/bwdmqZ17uV1DgP+KYToW9dS+vUU9T2vpmkhIYRvL8s9LDCmbgaXRqBECNH3OQ8Ftuvvdxc6tIn+w8uS/+I6/ga8DpRompYJPA6I/6I8gyOHemDoAIunjcjGtYehQBJo6bNtIH333VaP7NF7+rycmqbdo3+XLYTw7KGMvmXN2Kksm6Zp25F1qbf+CCEcyOmbIwajoR9cliPnx28TQlj0haZZyDnOPfEicIW+mOtADj/3FTeydxQVQkxErg0YGOwNK5AN5T1CCKcQwiaEOAU5BXmzEGK4EMIF3A38fYCe/+54HpglhDhLCKHoZU8VQhRrmtaEXHR9VAiRpdefU/XjWgBvz+KvzuPAr4UQwwCEELlCiHP1714CZgohJgshVOD/OMLaviPqZg80mqbFga8jF4/agUeBSzVN27AXx74N/AE5jN0ELNO/iu3DpVwL/J8QIoD8h/HiPpRhcASiaVoK2TkZCWwDGoALkWtIzwEfAFuAKNLY4MuUXQ+cC9yBXAeoB25lR7t0CXKefQNyzemH+nEbkP9oNuvWOkXAQ8hR6wJd5x8BJ+r7fw5chxzZNiEXahu+3JM4tBH64oTBQY4QYgywFrB+yV6TgYHBEY7Roz+IEUJ8Qwih6maS9wJvGI28gYHBl2VQGnohxNlCiGohxCYhxI8H4xxHCN9HDmlrkRYE13y1l2NgaNvgUGS/T93oDhY1wJnIebCVwEWapq3brycyMDjAGNo2OFQZjB79RGCTpmmb9cXIF5ALLgYGhzqGtg0OSQajoR9Cf6eIBnY4CBkYHMoY2jY4JBkMz9iBHHF2mR8SQlwNXK1/rByE6ziE8JCfYSeR1jCbBBb9ldY0tnU2AiAsOeTaLaQ0DX8shcUkiEa72DdryyMPTdP2h4OYoe0vjaHtwWZvtD0YDX0D/b04i5FedP3QNO1J4EkAIcSRbePpmMkJYwvwWs1U5jq5YHgWeXYLDcE4z230saChiykFbopdKstagixuDjA600YsrfF+1e++6qs/kjC0/WUxtH1QMBgN/UpglBBiONLV/9sYnpi7wc3VJw1FEXDLMQWM/N8TIKVBJEGpRWHs75cwMc9JayRBJJnGazOTYzVT2x2jLMOK/AkNi8sDhKHtL4Wh7YOF/d7Qa5qWFEJcj4xCpwBP655pBlRCwQRofkr/nMsH3/8xNf4oVx5TAHfMAnLhsfmsWd9KxU+mYjebSGkaqmLi1To/JS6V6cWZ1HZHaY8lOXfSzVhNAsUkmL/JR3mmDbdFQTEJ1nVGcJhNtG5eBqz6Cu/78MDQ9u4wtH0wc1B4xh4Zw1sFuIAzji9mtS+ML5rk3hNLsJtNTMx1UuJSKbpvBjK09w7euvJ+VrSGCCfTLGjo4vLyHMo9NlojSXzRJD9aUgdtPqAbaWpvQkZ3tet/VbBkQJYDl81MkcNCzYZPgQ8P5M1/5eynOfovjaFtQ9uDzVc1R28wIFbIc9MVSzEu20GJU6W2O8a4bDsr2kLE0xpFy+vhxP6VoTzTxipfmPpQkopsByUuFQCPquBRFSYWZbCirQUZSiekH6UhQ4To2QkTVmiFIGZqcAEZoFwEeU5oWgxUH5AnYHC4Ymj7YMfo0R8wzmHM2DI8qkKe3QJAXSDGlAI3q3xhnjt9BKUPXL/7IlZ8QN3fP2N5W4itgRgui4I/nuThta001W5BOtDGkMYhKWSvJ61/VpHznUnApm+LI9cWHWAthNhyZPy0ww+jRz+YGNr+KtkbbRsN/QGhlF/MvIwVrSFURVCeaaM1kqC2O0ZKf/5LP/39FxzbhEyWs4O6OQ+zuDnAm1v9lHts2M0mtgbiPPnRNgg3Ax3syDtiQQ59Q+yoJCqyQvjZUWFcQBEUlkPTn/fz/X/1GA39YGFo+6tmb7RtBDU7ADx78dXUdEUpcalcMDyLfLuFsVl2yj02Ch0qSz/9idzxj3/d5djvT9s5DSeUXlGJLyqtEe76pJGfr9zO0pYgE0fl4B0+EpkDIpfpJ5zGTdOngXciUvwmIAPIQmYeNCGTTVmQFaIemj6GwquQkZUNDHaPoe1DA6NHP4iMH3cjPzm2iNZIAodiotAp7YeTGpS6Vc759ng4c/qOA+6bB7dd3vtRCAEcjxb7G6ijdip9K/ll99IajkNzM1LcKaTo42Ap4NfnjMZtUVjbEeHJRTXIHlRPhegGuvTj0vpxCf17M1AAuUOh7elBeDIHHqNHv38xtH3wYPTov0I+vu73TBuSwZKmAKm0Rm13DI9qRhGCyhwHDrOpf0UAxK//o1cAYNG/9a3bmXbK4wOcYRgttY+iPTSLj66ZRt6I0chhawwwo+S6SKUhntKoC8RAsSBFHkZWAgXZ+8lG9nwEskKoyErVCm0d4L1i/z4Yg0MeQ9uHHkaPfj/z3VN/zPMb26GpFXK8zBzpZVK+i7EeOylNQxGCEpdKZXkO3PjdfseKnO+B7y9oWkJ+zroK/M8A+WjRxWDduefTH2H9DsSDQBeoQ7h6cileq5kXN3dQW++HcCeyZ+MCNBAqaGEgqL8UZKXQ9M9uYAiVFUOoWvPQ/nxMBxyjR//fY2j74MRYjD2ATD72Zmq7YzTVfoIcOhYy/YQJlLqtqCZBayTB5EI3UwvdFDpUcu7ZNbS8EGcD76JpjUAh/OUfiO99Gzn8/BaatlMGwM1ViLKfoGkL5OdQNcJ1jX7+bMgoxJ5poyzDxgSvg+e+PgbmXNq/jI7P4ZlP+FvVdr6zcBM0NyJ7RWlkxTEDxZxx/DAWfXzouqQbDf2+Y2j74MZo6A8Q00+Yw4KVm5BpMxuRw8MyykYPY5hbZdFGH2UFbiYXuPjVCcWU3H8xcmjZHyGOBVahvfH/YOY5QAAhjkHaEcfRtiyF0so++/f9fbPRNB8QYNak/6Miy8Hdj3693/57xSdLEBMfhVTP3GgSsELmSE4b6T1k448YDf2+YWj74MeYox90TOC5jAUrP0S6YW/Wt1uhMB+vzcyimnboilLsUvHFkpR8q4KBKgKxjfI4HGDq+d3clB81E8gF4ojhP4a/viK/evn1nQroQIixgJs3lt3P3W/98osrQtNq3Qoitet3x01Gu3s6ZA5FLmzlAC7o2sL7a1uYOP6He/lsDA5tDG0fThgN/T7zDeAC8H+GHE5akaItlq+mZlasXg3d2yDLzsRcJ28smwsTTx24uNYgPW7d9/1pZe/m6vUPAOX6p38jvvso3P8MdUu3AWfvVMh6hBA7vTIJ/PRxhJi1Y1vR1xA3/hI+WTbwtdx2OdqqG3jhkjOorBgNzjzABbFtVLWFgNP35YEZHDIY2j7cMBr6faKYwrJCoA5ZEaJAEb2eeCSRqV6DgI3xQzK4776zAac8fHPVrkVuaEM6eQzl7k+b+nyhMmnCmD6fP2X00x8z/G+rOHdSxV5cazcZd18DvNlnWzuQDcdN/uLDSiu5cGIJL505khlH5yOl4ibV2MYLl5yFET3jcMXQ9uGobaOh/9I4mX3KxTTVrkIKXiB7Oi49LUVYf5kBJ2rxUO49sQQqTgbk3KMoO36nOUigJYjsOZXx8flj+3219Nlv6e9cgJdSt5Vrjs5n1lAPT8y+a5/uQtP+vuPD1k9g6Xvw6N/673T9xZRm2PjDKUOZfsJo/V5D3PjhVv7wzV/s03kNDmYMbR+u2jYWY78UQ5k04XyWrVqJdMWOI21zh0LuEFSrmXh7CKKt+v6ZYLWibb0N8o9BiFOBxb2l3TT9Th5891fyw5qliGPuBHxoG/4Eoyf2P/XcPyGeWM74HCerW4L85xtjmfrEtYAdIUYjc1bvHdojf4VrZRh1Ic4DYnx83VnEUhoTvA4cd32/z95buXzq48wuy2Z9Z5QfvbkeEm0UlpXzw4p8bn/10KgUxmLsnjC0fThr2+jRfwluPfsylq3agqwIMaRtbgrwQDjB5eU5jB/phcJhchsmtlx3EuQfo5ewuF95Dy3o02PZ1IF08MiBAveuJ587De2hWTwxpZRHpo5g6mkjIFCHTGH6JduSnoqgXARsBTZz/KMfUdUeoqYrCvNe6rPzMB6dXMpdnzRiN5u46fQyQKGptoH2qJEU4nDB0PbhrW2jod9bcq7g/nc2ILPJ9djgupCOF0HGlWZRlmGlMtfJuaVZfHbDKWja3/tE7Uvvvvx8F8Uj85l87Hjoig6wwwgY5cUXS3Lda+tY89E26I5BdRPafT/9EjfSJ85HegNyDlYBbQMPrWmhuitK3ZpmYj97snc3x13fZ9nmDq57bR0PrWpC9vSi3P+vTbx++X1f4twGByWGtg97bRsN/V5wWuUtFHvswHqkw4UJKYggPe7Vl5fnkGe3cMHwLF598htU/OEmeODZPhYCyoBl985nrm6iNZpk8Q2TYOixA1+IQ+WccQWoeS6Oee5TKHTLYfCt36T5tkf26l407Tn55h+vI3tucWTljlBbvZqnN7SxpiPCe03d+O94DB54Ft5+hz9MHwV2C7RuQPb04pBq4toldTx78d17dW6Dgw9D20eGto2Gfk/kXYlDMdGwaSNSOGb9lUJaIVjwDs+lMtfBxFwn5zz1Ixg3ibo5DyN+dNnen+f8o4k31Mucml/ApieWI255mV8dP4S8bAeE4vo3TvLvvZY/X/jrvTiRFwBxzavIitwTuzsNtLJgZQ0LG7vxRZPUBeLUNwdgRQNlGVbunDwMrKVIZ5MwEKBhUyMPf97K9BPm7P29GhwcGNo+YrRtNPR74Nbjinh7xafIXNAm5JDWCmTqf214VDOBeJqxD30DKa52hv/uhr0q/9fnzpVv8o9Be/QSOHMkO7Lp9KfULcOu/uzj7bR2RuAv/U3Zrnzhf7nlrDv3fNLOz8FXo9+PgqzcHiAfSPLQR9tY1hKkJZKgPhin2h+hyKFSlmFj5rGFSJtqBRkCto0Vq+uxmw0pHWoY2j5ytH3o38Egcuf//Jz731kONCMfVRYyUcIoYAgoIwAHtdUtzBqbh7Q1Bu55Y7flaprW+7qj78r+cUXgdQC+PnvHePCCX0FgPW9v66LzJ2fw8TfGErr2JMjoCb/aQy4PvPOrPd6X9Zg/IBfdFGSFNgNp8JaRN6KUYQVunqpu47F1rdSH4tR2R7GbTSgC5hxTwMyThoJ1BHKIrwJpXlu2bkfFNjjoMbR9ZGnbaOi/ALX4at7e1gW0IkVjp7DsWK6ZVsH4cQWAH1JbgQCoFt7b1L7jYNO+WPJthm1d4HJC1VZo+QyCG/CU/pCbX3kckXElNy/bhqc0i4rTy3DcdTacNxY6BzI9m7LbM8UbupDWDHHkgpUZUCjOsvP6WaPY9O1jmFs5BK/NjC+axGoy0RpJUORQKXVZuajMy+zjh0DuSKQDTSvQxU//XcutZ/9sH+7d4EBiaPvI0/YeG3ohxNNCiFYhxNo+27KFEP8SQmzU/2bp24UQ4g9CiE1CiM+EEMcN5sUPHlP57igvVWs+pydmRuawSu48rojHFlazeu37yIUrB5CPd0gGP1pWLxd36OKBj+p3W/q0geb8Hv0INA2iEepf+AxRcAfC/U26tj6ODCa1jNrqpTz47kZinzTCwx/Cgo2wbBt8/hHQgvQ+BE374AvOfIr+tydMazbS67EQso5mekkmJ47Jw/yDE7ntnz/nqZ+ezjklmXisMllzOJkmkEhR4lKZMdTDvacMA/KQXpFZEGrklS2dXP+1L2Mp8dVhaNvQ9uGq7Z3Zmx79PHYNPPFjYKGmaaOAhfpnkPZNo/TX1cBj++cyDyy/P/8Mnl5ejxRXAmwVPDp5GNe9+jlSmHa5nS4giC+apKqpG+uDHwLdzHll9//5+4dEDcF98wg3BuCoXHh/Cw981gyWDL38vnzCvJp2WqMJkm1BuSmWkhXirU/gly/Br2XWnMKya3c5r3bPVRBYDzYP3z11GrecdTqVFRMYM3YslcWZPHV2OeQ5dxww7UxKc53YFRPxtIZbNeGxmgkkUpRlWMmzWxhWXgCmbPJGFENeKbXVTbgtA1thHITMw9C2oe3DU9v92GNDr8l/oR07bT4XeEZ//wxwXp/tz2qSjwCPEKKQQ4wVbSEIb0WaaI3mtzNG851FtZBoQg4H24EAUqxt0LQd2muJN3bzwPnz9FIu2MNZdHviTz/lyndquHZJHbfd/g7T7lxAayTBjGOHMGPit+mxJADA8V3OK83i5c0dMrNOdwyGZsI3x8E5JxBOpLj8XxsBaNz0CNrGFWi+tXLO9HfPwu2XM/KERykuzmRqoZvzSrO487gi7j6hmJsq8mnxhSEYh7rOHee8fhJ5dguqSeC1yuFuvt1CXLegmDXUAxaF1s2t0CoTN//mP7V7t3D2FWNo29D24artndnXOfp8TdOaAPS/efr2IUh3th4a9G27IIS4WgjxsRDi4328hkHhjONvYf7iT4BOIIebplfwo3droKkR2IZcTEoih71Bfb86oAvSAe5Y2UD7jx9D0176gjP00AZViznxipcJxlM4zCbe3OZnWUuQukCMcCrN5EIX10z7PpUVN+Ed/gNmjCtgcXOAQCJNMJEmsKUDKqeAZyyQQ10gxjOr+wSNGnkCZB8NBODmS4AktdWbiKSkg4s/lpJDV6tCWYaNpnCC8JZOcFt3lOE6ipwM+dltUYil0sT0iuC2mJjgdciMQqSQPcLtEG1mfWf0UB3mGto2tH3YaXt/L8YOtFIzoPGspmlPapp2vKZpx+/na/gvOI0ih4rs5KmQU8pTG9pAVZD1usfO2ITsEaWRFSOtv7rwqAq5T8v6XTb6+gHO0UMJF970Jo2hOC26u/V5w7KYmOeiLhDHahL4YylSGnhtZsLJNCUulYg+l7hwezcvb+lk+bW/h/UfATD2oRt5YdaYAc4l3c6FmA2E8G3ZQn0ojqoIAvE07zUGmL/JxwOfNTOvph3G7RQ5cJgHVTERTqbxx1MoAlRF4LIoBBIpJhe6mXzscGQALDPg5+0V1frB3+IwwdC2oW0OVW3va0Pf0jNs1f/2RDpqoNcOC5BPp3HfL+/Acuf/nMbzn2xHCttNYaadIocFLApS/Cakna1Tf99Tz03I4WqE1o4w+CMAbNrw892eL5xK09Aewm0xEUtrlLpVphS4eOnMkZRn2nh+o4/GcJwFmzuIJNOs6QizpiPCen+EQCJFOJlmvT/K3379HvzuWQAufPZ2WXjn5/DKG70VZc7ZP0OaixWAI5f3mgKs7YjQGklQ0xXl4dVN1HZHeXJ9K3y6qv+FThneWwnNQhBJpQkn01gVgcNswms163bQPa7zXUCYhxfVcuvZR+3rz/FVYWjb0PZhp+19behfB3pc4y4DXuuz/VLdQuEkoKtnGHzwMw1FCAgGkYtROTSF46iKCdr05AumCWAdibRI0NjhlCGQi1gWXBk2Jo8rgGdfRlacgZgAQI0/yuyj8/GoZrxWM5ML3Mwc6uHlLR14rAoTvA7Wd0Yp9NiZXODGH08R7Aj39j7iKSnQJc0BnlrRAFV9AktlHQ3nz4IxJwHQGklCoZuJ44dwUeUQRmfaWNEaZM5H9by1rYvvjctnTJad1kiSn/3sX/0vV3dbj6TSxNJpfNEkDrOJWErDZVHIt5up7Y4irSNakJVuO6SD3P/OOn2Ye8hUCkPbhrYPO23vjXnlfGAZMFoI0SCEuBK4BzhTCLEROFP/DPAWMufYJuBPwK7L4wcpt5w1ibfr9QUoFHA5IJnGahLI+cpSfnvuWCaP1TPS9GaU73mEGqAQ7IoyMc/J8o+2AfCLmbuGOm27/WoAbqoooC4oXb1ruqK81xRgfm0HJU4Vf0yavk0vzuCHFflUZNspcaoU5rtpCidQhPQoV4SgPhTnqg+2wL82wYuvsev6Ijx393TuPbGE2SOymVzg1oeoJibkODCb4M2tfmq7Y5xX6uGu97bA86/sUoaixy7p+dsaSRBPa/owV/YUZY8wjuwl1gFJHn5/M9d/7QKkh+LBg6FtQ9s9HG7a3pm9sbq5SNO0Qk3TLJqmFWua9pSmaT5N06ZpmjZK/9uh76tpmnadpmllmqZVaJp2UC1G7Y7KHCcr1vbY63rAbkGxKKzpjAAO/n3VqVx/dD4XjfSyI1a3CtiQw2E7oEGXnyXNQRnWI7mJuW/M7Xce19Dvk3PPbACuffGnrNjSQWM4ToNeKcoyrKiKiRKX2jtf2R5NoghBocPC5eU5zB6RTb7dgtfWJxNOMo34yd9g9rlIG+KdOOV0yjKseG1mIsk0pW4rZRlWJuY6Ge91YDebCMZTrPKFIdTKyQ/0DzurmkSvaVksrRFPaXhtZkqcKilNY60vBIX54DoamYvTjIwZ0gCJEA//ZzPfO/2aff+BBgFD24a24fDU9s4YnrE6td0xUHoCIWVAWiMVjhNvDfLCJROZVjkExQTxVM/iVC6y9xPRP0fpCaK0oinAitYgC6+Ro/6+buGBrY/T16xM63qGeEojz27BLCCSTOO2mLCbTVTmOlF0T0SvzUxZhg1fTC5uJTUNu9lEStMoz7TxvXH5QBZCzJaR+e75C9w3j79fei/EpVnaBc98hzy7GbfFhMuiUOq2Uuq2clZxJmM9dpKaxrINbRSPLGPmMA+y1yJxmE0EEil0owYaw3EUIXBZpIQqc10QTVJc0NPzcQBl+vPZDKkITy+v59xJt+7fH85gjxjaNrRtNPQAXMiyliDEw8iejHwshV4ns08o5rxhWdASZEVriEJnT2/Hi+zxWPQywshhpQnCcVa0hfDHk1Tf/Mf+p/psKWxaiUyKABBh8dxp/KKyiFhaI5bSeK8pQCCRkgtBLiseq4LbomBVBMPcspdR4lRpjSRk9EG9x1R+VD7kunmvMUAslKCuJcjq9jCLb+yJT1LMOUOzyLSaseq2w2UZVuxmEx6rwsn5LoYN9VDuseFQTMhAUJJYWiOlSceSrlgSVTFhFoKUJofZ4VQaEim98mYhh7cB5PqlAJoh3CkTMOdcsf9+OoM9YGjb0LbR0APw3VOH8/ZmH71R/CwO3ewMrh6Th1U10RiK4Y8nMQuQP3YIORQW7EgmHAN8EAwyf2M7VW1hFm7vll+Fqvnjt+7ivB+8Ch/VA8P0Y+zw9VP4ZVUj8zf5aI0kqGoLEYynqci2k9I0HIqJ+lAcl0WhxKn2OncEEinsZhPFLpWxHjsTvHLutT4Up6ZLWkeUe2y0RhLUzXkY2Aq/uIoLL6/EazMTT6dJaRon57uYWugmlta4uMxLRbZD713Ze5/R6NIsWiJJFm7vRhGCrYEYLosJXzSJP5ZkWlEGBIOYhQCvGznEjSIDS+UhE01HaGjsZvrwbA7HBMwHI4a2DW2D0dADSBG1NyKHqXmUjfBS6FDJs5sZl22nuiNCayRJIJFGEYJxRxcjf2AbUMCO+UwVOZxrg+YOfvNONfetbpLxP5wObpj/bWaPyIbvnr/j5J8sgbpaVEUwc5iHlKYxNsvOOn+EqvYQpW4r6/1RSpwqwUQKVTHhtiiEU2kiuhlYvt1CsUul0Kly2ohsmkIy6l8gkcKtDz+Xt4V4+bIXuG3Gz3nrr58STqYpy7BhNcnvFSHwqArL24LU+KP6MLaP9eBN32WC18G4bAfFLpUy3dHEH5c9oNZIAohTW90MkQTYhoE6FnK84ChEDnW3Qng7C1bWc8bxNw7Sr2nQF0PbhrbBaOgBp/7D6tloTNIBozLXQWWuk7UdEVa0ytgbDsWEqgguL8/BXuIFhnLZaSciHSTzQYwH73HI//RtoLWxdZufRR/XM+2EB2Grn4uf/3H/0x83GUoree3jBnzRJCkNJuW7mJjnxG1RqMxx4I8liaXTjPHY8agt85S7AAAgAElEQVRK74KR26JQ3RWl0GEhx2rGH0tSF4hRH4rrXohx8uxy+O1RFVwWhZEZVla0hgjqtsoeq0JdIIZiEsRTctgcS6eZXpIBr3/a71JHjy9kdKaNfLsFl0XhyQ1t+KJJmsI7HGMgQXFRBtjMuArcIHQTCtUNjEMOeZtYVLV9EH5Lg/4Y2ja0LTEaek6iqr1nqBpHLcogqWnUB+NMLXRjN5socqhUtYVIaRqKEBS7VP4+bSTaS9cw77ZT6fEofOJbFTx71ihmn1KJNFvrgmgr5LpxWxTWPPQh/P3VPudOwzvv8MYV90NCxgc5rdBNPCXFntQ0WiJJhrmtmIWg3GMjmEhR5LDgtig0hhNYTfL62mNJKrIdNATj+KJJ2mNJWiMJ3BaFlkgSt0UhpUknmJSmUR+KE0+n8cdSvfG43RaFumCc0Zk2wsk0yZ0Fe/k3GZdtpyWS4OkNbTj0gFCN4QS1XVGkGV6ShlCczEw7Fdl2Mh0q2MyyUjjcyKQWEdAS7Ig4aDA4GNo2tC0xGvqCUn0opwJllGfa8KgKRQ6VsVl26kNxYuk0ZZkyol0wkcKhmJiU74ILvg7nzOCy08YwY+J4xmTZmVqUwR9O6ZmjVAEnf50+irmVQ6SH3YU9MbIS8MunYWUDXpuZa6aNYXSmTR/CSmuEfJsFh9lERbad2/75c0Z77DSGE4ST0jxgglcONf3xJPl2My2RBKlgjLfXNvPTD7eyoi2ki1zONxY5VKYXZ+K1yjnEedXtrPKFWeUL815TgOquqIxW2B5mWYs+J9n5eb/HZf7l1SxrCTImS86xuiwKYz12rIoJlExwZkMkQU9SnopsO6TSIMCb70I2HG6giZknTRq0n9UAQ9uGtns5+FYNDjDDMmz6O+kFKEVjwW0xUZZhJZBIoZpkr8BsEtjNJiLJdO+PDbCgoZvbxhfw9rYuAokUhU4L8tHmUH5UMVZFEEulcV+n//hz/4Q/nqQxlKApnKC6K8rSlmCvs0ZZhpUSp4piEqQ0jUuO0r0Q7WZq/FHCqbQ+zE3jtih6BD7ZA7ll0lBSmnT4+HttB8PcKlMK3NT4oygmiKXMzBjq4cn1reTZLcTSadZ3RnuHuPXBOF2dYVwWE1XtISqfWAE/PrrfM5v7xjX8/dJ5rOmIcJrNTFcsSXs0CV4HY3IcrO+M4GsKsCyaZFy2A4JhAHwRGzJwlgNQWdDQDWQgF7MM9jeGtg1t93DE9+grsu3SfIokIFfre5w1PHYLqTQ0haWJVzKtkW+3UOhQ8dgtvWU0bnqEH35dBlyaWuTmp8vrASfYMnBbFBQhKMuwQdF4Fv/gd7xV76eqLYwvJoehL27uIJhI82pdJ4VOi27WJecUbztnNKXP6/E5zi7HH09R45e9E180idcqHTtquqI8v7EdfzzFQ2uaaQwn6Nq6lQUN3TSGE4zLli7gdcE4b2/zk2+3kNI0HlrTwvPrW6kLxHi/uo2u9hBer4OKbIfsHSXTSPO6vuRz4bO3c/HIbHzRJA+ubcFhNlHotrK+Y8e+LovC2tYgFGbJeV8hkAt6UcBEPBADTt//P6oBYGjb0PYOjviGvtRl1Yd7bjDZiKXkXOXEPBf13TFaIglcutfc2Cw7Sd3mdlNnpH9Bl32Tu9+6ihdrO6B5G+DBleeiMtdJTVeUnLvOhOQmPKqM8+GySCeNUpfK+/Vd1NZ1opiE3hMAu1manfG/s9la8wg8/gIcczK/nzRUDhmBYCKNw2yS8bsBu2Kixh/lpooC3v+0EWjn/arlxFMyfkhFtp14Kk2e3UJS03hzq58Sp8qdE4v5YUUBs8cXMq40C19riIcWbKC6K8q6thDJXzy/4z43ruh9e8Ezt1PiVPnhuHwqcx1MynfJXmRXlJkTChmbZeeio3Ipy7DJZ6xpyKw/NiAKXUFg9WD8rAYY2ja0vYMjvqHPs1uka7Q+ixVMpPDFkpRnShtdr82MVZHOE4FEiogebMnXuxLflxLdrTwNVhfBhi6WtgRljOvX11N969uUZdjwx1JYFROKENJrESDmxywEpS6VcDLdO3/Z480orrkIgJwhGdSH4kSSaVoiid5Ie/FUGqtioqYryvrOCKR9yB6GnxWtcrGtMZSg1C2dVEDGA1nlC+OLJlnTEUZVTDjMJjnvSAxfNEkkle5NxADAqIl6uZJZxxSwqTtGnt3C6EybfC76lECeXWbsqQ/FadjUBC1BKByGXLRy6CW07J8f0mAXDG0b2u7hiG/oS90qb9f7ATtk2VEVE7FUGlURNIYSmIWMaw3gtVp6ExPk9xne9uXeiSUgcigszuTqqSOozHFQF4yx5r3NRJJpHHqca0Wgm5xpeuhXH7WbfXrMcPldz7ymFllJj5fipi0d5NkthFNpyjKsNIUTvYtYLZEErb4wC1bWIZMlOIEED/97KU9taMcXSxJOpsm3W7CaTMRTaVKNjTy2sIZnPm+hqj0kky2MygGXh5e3dLKuM4LDqgDNfe5yh7PJ0pp24mmNsgwr47LtTC1yo+a7cVtkiNdAIk28YQtYHUw8ppDJBW5wZCIrBIwZe9X++BkNBsDQtqHtHo74hr4+FCdS7wOawNeBqsff8KgKDrMJxQRuVT4mqyIdL1Jp2QNi0b93KW/WX24CrYMpBS7GeGycnO+iQg+RWuxUCSdSeHTLAH88RW13TNrkEoZEO0tbg9QFY7gsMmEx/3gdbOVo6z6A2EZ+s6qJCV4HpS4rKU0Og3tEN7dyCHRFdTvoXM4/+TRgBLCdV5Yu5Ir5r3H2U29z/CPv8aN/fkDNhs9AuMBmhfbtrF/XSF0gRlVNO9jMlLpVyjNtVHdE4PH3Bnp6KAKuPzqP2SOyGe2xMaPEw03j8rEqJlIajPXYZE8nFmHF6nqWNAeYVJ4DTicUZkknG4NBwdC2oe0ejvCGfoQ+vAwiF6yaCUaTlLplb8IfT6Lq3nVuiwkUGeVOMUkhU+8foEwVLDksbQny4NoWFjcHmZjnIqlp5LithJPS669nCKsIwcn5LqRTi+D9TT7qAjHW+yPUBWIsXLARnnkJNrQz95t/ZUF9F/5YqtczMJhI4Y/L7PUTvA7OnTSUSGuI4pFezAJumn48ZJ2MDFSlAXZQHGDNBzpAq4Po50gX924WVG3HnufkpuNklrzVvjCKgHXrWyG2ccdtxjdSN+c1xmbZaQonaI0kSaY1fNEk6/0RqruiNIbjNIYTMohUXhZgJlM1y+mESII8u4U3tw70DA3+ewxtG9rewRFuXhnX5+ikSRTIOTig14EkkpRD3UAiTXso0bt94HlMSXmZl5paHyR8bG0pZNqQjN5YHVZF9JqzeVSFuN2sl2UCEtD1OR71DOxmE2s7IpANT71dgz+eYllLkCKnyoq2ICVOlRKnisuiEE9rhJNp6gIxKrLt5J8yjDUdYV6saSfPbWVyaRZL0hr2jKHMGuphmB7V7/blOVStqUHGKU8CHop7wsRazXisCuv8USblu4inNb5/zjzKM21MLnDRoMcnaWpK4FEV3TFH4lHN+GMxFEX0up+7bGasw7Moz7RRH4zjGZNHIJGi3GOjatB+3yMZQ9uGtndwhDf0DXpoVh8yFogCnWGKHBbWdUbx2syUugWKSS7k+OOyApS6rTL+RZ5rwFI/PHcsuffOB0zYM6zc9UkjT0wpJZxIUR+UrttV7SFiaY0xuuOKXJhKAjZe2dLB5eU5VLWH8MWS5OnmYuO9DhxmE1OL3ATiadyqiQX13dSH4nhtZm4bXyA9sk1CivO4IlojSa5dUgddHUS6LFx0RlmvF+EvKoewsNDNwu3drN3SQWaui6lFGTz3/m923MwHi5hz93+YXpxJZY6Dv2/uwGuT4WCr2kOsbg9TmStd2nsSMhc6LNSH4gQTaXKsMjzt1CI36zujjMmyU5HtoNip8tDaZuZ/1gxMQub/MNh/GNo2tL2DI7yhlyFKpWODnkWmwEVtd4yyDBtui3Qg8ahmFCFzSEaSafyxJFMKM2DG2bsWWL2CnInFSFOrbiL1XXzv3DHUh+IUBmK9K/X5dguKSPLk+jZaN7cgh7dRwANNDQQSZYzJsmMWUOhQiafSbA30XJfC1Ce+BQxlKsicmtEkpNKE9aEvQJFDJly+aKSXp+tldqFYWiPfbu4dIlfmOKjMcbCqOAOH2SR7ZzQBhfJ+Tj2DB9qC/PGFz4inNKwmE8tagrgsJprCMsDUe40BCp0WyjKstEQSveXU+KMkNWntUeq20hRKUOqysqQ5wJOL3gfPSGmG5hoJwVp2pGc12B8Y2ja03cMR39BLTEAYCqQAfNEkz21s55wSD5PyXQQSqd6wpz3DSVxqvxJabn+UgvsWARquobkUj/SSby/oDR61pCnIus4IJ+Y5dW9EwfxNPuZv8iG95xQgAGIE5Lu469NG/jl9FP5YkhML3VgVQUW2g9uX13P3W7/sf/m3XNr71qG/ztE/TwAuAJ7SPz/17bspclgoz7RRF5Tmb02hBMFEutfKYfLt/yT/3j6Z8i74Or55n/RWsvqg7GWFk+neCqEqAqtJmuq5LQrLWoLY9YQOFdl2PKqZcCrNKl+YBSvfA6rAH4bcY6EtCJwKvLTPv6DBF2Fo29C20dBT0xVFCtGLSzXjMJsIp9J4rWaWtwUpdFo4Od+FP5bSHU7S0iU61n8es+ChpeSNyKO1I0xwWx1BPDSQZtMwmXGnqzPMApeVK4/KRRGwyhdm2apG5BxiC3J4K0DrojK3gD+fNhyQ2W/q/FH8sRTr/BHeXfHbfbjLADIGB1z5wh1y07Mvs+T/behdOJtc4KKqPcyUAhev1XVy9a/+DD+7qs/x0omlpxKEk+ne1G31oTjelBlVkVYSVkUwMddJYzhBJJmm2KViFoIpBW6eXLQRemcuP2Vc3hTWdsmgV8RHIlOyGuwPDG0b2u7hiG/o12/zAwpkOvRUZ7rpF0lUk2CVL0xZhpVUWpp72c0mzirx4A/GetMBO4Z+n9PG5TPB60A15XD/O1vAko23OAtfUwCinUCSVDc82ViH7OWYkOna4siMPibAQ9noEsoyrCgC8u1mqv1RKrJl3stpf5qzD3dYr5/HJz/WboGyoXDpBVx92ieceO5zMtuPW8bn8KgyAcRbdZ2smDWXuVdUUr24jhWtIcoyrMT0lHA9IWAXrrwdGb+8h828e9UrtEeTvXO/bouC1dSTmu0f/a5uYp6TukCMYEcY4qMwGvr9h6FtQ9s9HOHmlUBwM3ic0BXFZVFojSQYrUf581jNVGRJB4pSt4oi5JycP5bAc0ZZbxGR1hAXDM8ikEhx/zurADdjRuUwc6iH047ORzphJJBed0HIGctFU04lb8RIUIcBRUw+9jjyRhTREklw/bj8Xs/ClAY41f7DzS9FCTJrkFe+yo4H8mDbpzAsi+Wr5jK7LJs1HRGuPzqPCV7pBFPjl16IpT95l6P+uIza7ihV7SHiaY1VvjDzF9eycOUD9K8IACM4688/Ys4xBfjjSawm6aEYTqaZWpixy9X5Yykm5jlxZTvAW7CP92gwIIa2DW3r7LGhF0KUCCH+I4RYL4T4XAhxk749WwjxLyHERv1vlr5dCCH+IITYJIT4TAhx3GDfxH/HalDNQIwcm7l3Vd2jytX3QCJNayTJ4uYgwYR0EQ8k0uDrEwwp1s5Da1pY3d6zzYxDMeGLJaVX3dFDIKMEuSimMn14NlePyeVvZ5Tx+sUTeOnSY1lS10lrOM7UIjeKEKiKoLI8h5MfuxluvexL3lPnnncZeiz89VMgzuVnjmRJc4Cb31jP/E0+fFF5v3l2C+eUZDKpIh+72cSyNS3Mr26nMteBd3jObovPv/daphdnYtdDIaY0ei07+tKzLZ7WqCzatbIMJoa2DW0frtremb3p0SeBOZqmjQFOAq4TQowFfgws1DRtFLBQ/wwwAxilv64GHtvvV71fCVLolqnDFCHzUPpj0nognEzzWp1M+tsUjrO0JUhjSLpk13/cP3FBbXWLjAVSkMeYsUMpy7Cimky9YVnLizKYfsJJkDNCN+GSoVgVk4yERyIFKY2mUIKTR3oZ+bsb4Mbv6qXXw6oP9/J+fPCPxXu363fOB3K5ZN4n0iXdZeXFlQ20RJKs7QhTH4pTnmljUr6La8fmMfukEjDJnooeNny33PCPOeTbzTQE4ygCXtzcscs+XbEUOVZzr9fmAcbQtqHtw1Xb/djjHL2maU1ImyQ0TQsIIdYj84udC9ICCngGeA+4Xd/+rKZpGvCREMIjhCjUyzkoaaptAWQ41MZwgnNKMqkPxXFbFDyqmUZ9xb7EpZLq0mgIxXljq58eyWnaOwjXpcS7ozxyVjlFDgvhZJoFDV1UZDm449Vf6HuGuOHM3zBjaCZ2xUS+3YI/lsJlMTGxzMvyi8fDbZcDHbDo36z6p0yMMOGi8XDyVGTS5ggyOfEX8NRCcKoQ3ACuo/bq/lsjCRSToDjXydZgjFeWbgIS1KCydnQhM4d5iKU1LhrpZZjbymPrWpk1zANNq6Fw/G5KtvPbk4by4Npm4mmNNz96cZc9qtZshYphxNMaVZ8f2CBQhrYNbR+u2t6ZLzVHL4QoBY4FlgP5PQLX//ZMaA1BrpL00KBv27msq4UQHwshPv7yl72/+RgKM6AlSJFDRvzzqPJ/oMNswmNVsCqCukCM95oCtEQSu5SgBZ/lJ6eUMr04g/JMG01hmeps1l9u7bOXk+nFmUzwOkhpGnbFpCdBCPP6WaNkRXjoedbc+BwPPLwMu9nEhDynXhHk8Xy45otvo+NzyHXKyvDmBnjrbXjljT3efbFLZU5FATccnc/5J5SAsCA7u+3UVn/GQws2sHB7tx4JEaYNycAfS/G3W9/eY9mO7x2PQzHhj6WQ86k78wFVzQFp45yODPD9gcHQtqFtODy1DV+ioRdCuICXgR9qmra7tCkDjVO0XTZo2pOaph2vadrxe3sNg8dW+UcI6kNxWiJJ6kNx1nZEeuNht0Rkxpz3q+pZ0NBNiUvlj9+6q18p47LtdMWls8b04gwZBwMAfVj3++dojST00KoQSUlLiGdq2uWC1KaVvPzJdvzxFHN+fSaj754OPzuv3zk4ZTx8sGjg2/hPLXREIJ6EaJLqf22CxsAe7356cSYT82Qi6SkFLsaMyQOLh17Xeep5bdk6Fm7vxqoISpwqgUSKqvYQ1Tf/cfeFt4eYUugmz27GNfSkAXbogLZW/LFd5zgPFIa2DW0frtruYa8aeiGEBVkR/qpp2iv65hYhRKH+fSE7XL8akMvhPRQDjfvncgeRrihoGpFYilK3itsiU6ylNI21HRGCCek1iFVl7edN2M0m1ndGePeq38KfX+TR2b+mLMNKezRJayRJPK3x/EYfJ064Gf65FO56ikeXbaMuGMNjVWgMx6kLxKjtjrH85e8A8MefvMsFc6Yw5fFbYMxJYB/NLkPZziYp9nfe6b99zVLY1iWbos4oyS0dXPX+FkY/sozwnU/Ah//5wlu/8J6z8UVlwuWKbAczh3qYPqEI8CDDwZqBNpZ8uoS7/rMZVRGUuFT88RSjx+bBvJcgUj1w4RNHEU6mWe+PMrXI/QVXIKjZ0IL0ntzdcHn/Y2jb0Pbhqu2+7I3VjUA6n63XNO13fb56HehZMr8MeK3P9kt1C4WTgK6DeQ6zl/AqsCjQHpIJCpDDPrdF6c2YM3OYh3EjvUCSeErjNH0lfWnVdlb7wtyxsoHbl9dTH4pz8mvrWeULMzHXyapFtby1pYMJXgfTizNpDCVoCMYJJNLMOv9o3SwMbjirvDfw1BeyqBa2d0vLiL/8A/7wPLz4Gqxuhmw7pDUIyYo2Mc9JTU++Tq/ji8ssquCxda2Ek2naY0lGZlgpy7CC1wVqLtJW2QlYIbyFqvYwXquZplAC/nciXF5J8p73v6DwQrxWM+Oy7f2TPPQjhbSFDiPb0dN2/wz2E4a2DW0frtremb3p0Z8CXAKcIYRYpb/OAe4BzhRCbATO1D8DvAVsRnoH/AnYVyPZA8xaMgvckE6wpiOCP57CapJOJEUOFdUkWNYS1F2lZYKFVFqTyRsSKZ5ctI5FHy9j9dp1XPePVeTpMa8DiRRvbPWjCIE/nuzt6bREklx+TAHMmrHjEpoDLHx8xYBXd+UZd3DHOb9g7rxP5GSBhp6nEgjFwR+VFaElCHqWnlnDPPx5+iiueH8zHFU2YLkS6XBy81vVvL3NT6Pu+u1yquQVZ0JeMbIHlA94WNYSpD4UZ1y2Xd8+HPOxRVRd/+CApZfpSap/cmzhF5z/LaSHYidybfSAJVQ2tG1om8NU2/0QmvZF/4kO4EUI8dVfBAAzoVA6NhQ7Vbw2MxO8DprCCdZ1RhibZSeQSNEYTjDWYyfPbmZFW4gxHju+WJL3qz6HzAJcmTYqsu2MybIzf5OPP586HMUkMAtY1xmlLhhjfWeE6cWZzP3aSCjPgRmVkOyCLR16SrP+jB5zAzUbmsHjZOMVlYzMtIMJGgMxitxWljYGeKq6DbdFIc9u7q3AJ+e7UEyCGn90QO9D/x2P4fnpGcz99nx++eYqcLiZfrR0alnY2M0Yj63XHbwxnMAfSzG1yC3vz2Pn4rvPknbLIHN/emzw7f5zr9U3/5EafxSPVeHUJ25Hzo0eWDRN+0ps3AxtG9oebPZG20d8CIT+rMFrK8bXFiLlUHs95eoCcTxWBY+qSBtiJYqij4VmlGSiCEFJWiV17FiWNAcocakUOlSmFrp5sbaDYCIlo/Sl0/hiSeZv8uEwS4sEr80MS+o4b1EtJfdfJy20B6B6/S0IcRPnjhlOMJHmvlWNtETkMDnPbsEfTxFPpfGoCk2hRG9M8bIMKxXZDplL85mX4LJv7ig0Uo3nhBJw9oSkNUMi1RsjRDUJPFYzZiHw2mSmn1RaI99uIZWWMU2afvgmc17RK8MPJhL72b+x3v9MP0eYFa1BvDazbp1wKvAFC24Gg4ih7SNZ20YIhH5spdRlBbdVD2kKVsVEkdOCWQhUxURZhkysvLgpyCpfmFW+cG/S4ikFbsZlO8izW2gKx5m/qYN7Tywh02qmyGnBazOzpClAxBdmRomHeDTBm1v9tEYSlOQ64eXXd3Ntw7n3vOP586nDebVO5rsEcKkyGbLbYiKWlqFWVT3hs9Ukw62u6QiTTGu8u7iOVTc8BGv1+Nj20fCNmbCill9WbQcEmEzUBWL44ymCeiYgVRGUZ9q4ZJSXKYVulrcGWdDQxZLmAC9v6WDs0TdSN+dh+KQRa4aVxrYg3DtPnmPFB6xoDVHjj0pPQe+wwfnpDPaAoe0jWdvG1M3OiAspLvOS0mSPoSfpcI7VjF23O17TEcGsz0uO9zpwWxTiKY1AIsWKthB2xURXdxQiCV74VgU5VpkxPpbW+PZzH0BmLt5sO3l2C3MqCvDazLJ3MioHbvjOF19bpBr+tJK5/9pEIJGiyKGyzi8rhVX3vvPHU/I6VYWUBnl2M2cVZ8pY4UHpDfh+U4DxXgdnVRTw4AdbeHlLJ0uaA9AcAC0JLpkxiKgfMDNx/HCmFWUwudBFXSDO0xvaUBXBssaAnEMNhYEEeLIoy3dT7FIpy7BS5LCwojVEfSiOQzH1JqVo3XzgHUqNqRsMbR/B2jZ69DujfcY5Qz1MzHP2ZquxmgTtsSSqInqHjj09gfWdURluNRgjz27Boyp0NXZBZ4SZEwpxKCYWbu/mG29VU+OPAhYKc5z4Qglmj8hmTJadmJ4ubfHnLbB5N8nH7Cp0RvGoCsGE9NN26OFTFSFTm7VEEhQ5LPiiSf5/e+ceHddV3/vPnjNz5i2NNHp4ZMmSLVuJHDs2OJgYnCYk3JTQhNAkBEJDU0hvumhouS1vbstabVlwuVACXKALL56BEhpKyMMkweA4ISHGTkxs/FCsWLZkyXqMNdKMZjSPM3Nm3z/2kWSCEzvEeo33Z61ZM+ehmbPP+e6f9uO3fz/LVjk8t5+YoKPaR1PAVOFYvW4e6B1n9Zbd/OTYOG1hLx9ZF1OeGWQhMw75FGq1oqAtZPLBtY289Q2tNPrdRH1uGvweNjWFYTJL+wUxvnLT6/jAJc2YhmDvaJbedIFsqUxHtY/BySJ79g+QKJRoCZovXT7N7KK1fd5qWxv6P6CLwaya2d+8JIxdlvjdakl3rlSmN1Mg7DGwnWTBz96wmh/++N0c+vBlDExatIW9tC6v5ZbL2rhxeS0NfjdXLa2iIxrgU88OAGW+vrmVb75pBXd01hP2uKjzuklaNvf3jsOPT10d+KIATl/dyUhWjVv63S7CHhf7x3NkijZZW8XObgqYJPIqRVvWmWQaySlvC5X4uUTSUmOViXyJpoCHG5fXsLrGDxE/4ENVgjwQprp1CTesqKXxNU0QC3NjRz2FssSWEltCrL2Ru9+0ArssydplLFuSKhRJFWwCbhcdER8rq7xAmcPJvJP1SDM/aG2fr9rWk7GnYe9olo0NIfrSBZpDJqmC8v9NF1XqtYJdxnAJDr3vEvjQX8LB38B738GusRxrv7OH1RE/t3XU4TUEXz4wwr09Y9iDverLPXX8+Y9+x4aOOjWJhPKCsCU8MZgmmy4w4xV86rLqBKWTGQJu13RrbP9YjrWOt4S6NoOkZWMagqwT42NqgsktBD0p1TI7lMzxSH+K+ECKh4o2mxpD2BL+8TUx2sJt/P1//5apzECpwRRbuuK889IWuLAeVjewekcPIY/B2hpByFRxubcPTrD1N0cJLasHIVhT66ct5CVrl50kyVnCHhcH+pOz/wA1L4nW9vmpbd2iPw0DRx5jdY0P01B5NAtlSUe18pctSclApsDvbloDH3L8hC9ylj/fvJYDB49x4/Ia2sJeOiN+7nnyEPbQVDKGAJsuaqR6SZisrSpWyGOQs8vEc0Ua/B4C62IqjsepPP5Ltv/P77L9xASJfImhbJGAobIFZe3yKTFLXFi2JFmw8RouvC5B2unO/rAnwdPxDDtH1Ct+ND71nLwAACAASURBVA5Wkdx4jvv7kir2SFRNtn3tHRvAaADcEA3w0XUxaI/C+tXQtoH2Kh8pp9LtOZnlBy8keHtrDQ0rGrmg2sfVzdVc0VTF+jqV8CKRL0FVo/LCmJjf4E7nO1rb56e2dYv+tDzPYLZILODBEPD4UBrLVunFmgIe+sZzuDsbgPrf+6sb3/UjjCaVSzNiGjQGTMAEmQR8EImoTD2GC7ssGS2USBVKRH1uHh9MY0vJ/T97nga/hzcEPJAvsevBLgYmLZIFm46Ij56JAn63i2xpppUDKjOO3+0i6nNj2WVMlyBrl1ldoyrxvT1jZAZSUC6iVum5AC9ksjzRFeeJ3nG+cVU7a2sDPD45QUNrhPhYlo+sb+JPWyNQF0C1hGrYUB/goeNJdsUz7Ds2zvUXLyEW9HBFLIxpuNhQF6AlaFKwJYPZItsOxdm0Kkp/xmIxRAyobLS2z0dt6xb9S/DtHfdQ53Nz2ZIwG+qCZO0yO0cybHvmBTY3VzO690Ur3/c/zX1P7+MLl7ZwOJVnd3ySkayFGg8cAzJsXl7LSK5Eo8/NUFZNfFV71f/arF1mJFfk/t5xlevzmRNs+c4e3C41EWVLSdTrVu5lZcmmxhBtYS+5UhlDCKJeFRs77HFhuAR+t4sLqn1sXhLmb3/dR+b4EJQTqAUdATCqUEsQbajysqmlmr95/BhXPPQ8//zsCeL9SW5Zu4Q7OuuhJgAH4/Dro0CGK9qjJPIlupJ5quuCeF2CA2M5RvIlEvkShpNaLeB28eXd/TAZZ+dAikShxHSQLc28obV9/mlbG/qXpJenRzIkLZumgIf+jIU1MAguFdPjnp4E3L91+mxxydeANNsHJ7DLkrDpoiuZQ3WavEBeTXRJSYNf+R3f+VQvdlliupR3Qc9EgUShxCPHk7z7oa7pNGWGUOK2pZok2zOaZXd8kuuWqVbUVDc35FHjom0hLzcur+GG5TXYUmLZZec6QoALXB6wyxAIEmtvYtOSMDtHMpApkCoUyYxOYjaGuaalmpXRIEzklavZWA5OHAFDdYU7Iz5SfQnu/fUxPrbtBQ6MZems8eF1CeWtkS5AcgLwwegEuYEU87UEXHMqWtvnm7b10M3L8MDOz9Pxln8madnsG8kAWSgLHugdZ1NjiEs+9UvWfOnX3N+bBGsSvMvZuneIBr+HpoCHTz/VB2RgyQUwPOaswLOJ54q0hUx270vwsV0m1znJD65aWoUtJelimf1jWRIF5UXQUe2jwe9hS9dJ7u8dJ10ss74uwN5Elp/cuh6AZH+S/WM51kcDhGNV8KYVUC6z8hdH+MDxFF+1JZzMgMcgtCRMW9jLgYPDtITUcniGTgBFyEbA68Uuq0q760QKr+HCHpygYyBFeDzHf+0boms8R89EgdCyWjK5Ite217KxXq1CNIRg/1iOPaOTqBRzKiE1VSFIzc+z1Pw+Wtvnl7a1oT8DbiFo8Lv5xCXNfHZ7Hoo2t11Qx1VOdL/eTIHbL6zjHrfLCQoF3/51H1hTEzMhPnhxjC+XIV206RrPEfYYbI6Fob6GgNvFYLZIrjSTszNXKjOQKcDwIPf0H0BF1wuiIuG5oaaaq9cv4ZPrm+BtnbC0kQhlLqOMypGhWln86FEu+8FeEoUSDSGTuOUD003ACUMLOXYPpLj5wqnx2AIwDgU/9mCRp4bTfOjiJTw9kqE3bfHzgRSZ54boiPi4vq0GryH44ZEE++IZLFviNQQhj0Es6GH7iQm++FQfqjstAR/RWj+JBVwZzje0ts8fbWtDfwY++8hv+cbNl7Lr5CTXvraJiKlW+vVmCrSFvaoFE/bymY3NfLd7lCf6U2BlIRCjs62GgOPdQKbA1r4kDX4PDX6V6ac6YBLxqhybXpegLeylYJf5ybFxGB5HjYF6gCSqW+gGApA0+eSmVri0BcaykOyDpVUQWT1z4Vsfw/uRR7BGJwEw64IYARO/20U8XQC3C39LjPYq5YGhvtt0ftMHLhf7x3JsPzFBf8aiN2MRMQ1KUtISNLHLkoLj3nbDa5fSEjJpDXvJWDa9aUuVYTKDavEod7YNdUG2HZu7Z6d5ebS2zx9ta0N/RrbyNw9U8+At6+hK5nh6JMNAxmL/WI4N9UGaQyZPDaW5ZWWU96yMkiuVydYHp7PahD0qWFRvukDIY/BA7zgBt0r8kOp7nl6jk1jAw3XLIiQKJT7z3BCZ40dQ/UA/KjlDCLBQUyolkMOIj/8Q8GE0LeEDaxr50p8sh497wWgHRrnkk9uwckXwuSE5gTWQh4CfcKyKL21aRtKyeXokw31dcQ4cHIFQBDIWqvIVwVk005u2prvku+MZNjaE2D+WxSpLIqbBta0RkgWb9ipVEUKmwWAyT193HHCGBCgC7mkvCs1CQWv7fNG2now9Gwr/ydu++1k6I37WRwOMFlQ6NoCoz83hVJ79Y7lp74H1UeWC1ejzEPYYDGRUfOsN9QFuv7CenrEsI7kSkCN+dJjulArvOpovOZlqkqguoem8T6Iqw9QEWL3zPoE9+AJf3rYT8flfwRefgtIR2LqbBr9HXXuuBIYPhBuyk3REfHzpwAj/cSjOfbv6YXwUvG4Ie53fE0ARTIO1tX6SVknlAHW7CHsM4rkiXcn8dHwTr0tF/4v63BTKkqFJi68eHAGGTylDNcSi7N536EU3NoaUEikl1L13Fh+g5iXR2j4vtK0N/VkzTrpoE/W5yVhlroiFyZbKSvghU/kDWzaGCxL5mbgXpktQktASVMu3OyM+NjdXk0pMorwE1ORQsqA8IFRkwQhqLNIDgRAQRrUecqiKkUNVjqJzzA8CxEcf4pNv+z4P//Qg3ak8l7fV0N5WA2ETgib+lnre3V5LxDSI+txc+ZoYm9avpLmlGmypfg8/iAauWd/E5iVh2qt8RL1uEvkS7VVeYgEPrWGTiOnGslUlaQubxHNFklaJTz83SPzo8849K6C6th7WRQPAtt+7o1LO+B3Lk9+encemOQu0titd2zp65SvkO7d8hj0nJ1kfDTCYLbKxIchAxlJeBAWb/knVHdzYECSRL3FgLKfyUAZNeiYKdFSrhSGPHB+H4TGIVPGBS5pZH1VxtUdyyl/3kf4kO/f2AwZ4PVAYRk1aeYARVNfXhxonnIRAM2ST6jMNUB9UWXlKZe7YsJRsqUxL0GRNrZ+dIxlsKTGEWnjyxGCanv4klMqYDSH+qqNOxSZ3qYTSU4kZ0kWbDXVBJ365aulETAOvIXhyKMOntnaBTKEq6ShT7nehZWvJHP8ZKuXqDC/WnhCbgV/PynPT0SvPjNZ25Wpbj9G/Qr7bPcpVTVWM5IvEc0X60gViAZVMuCPiw+92TVeMnokC2w6OQDYLsVrW1AamW06XL60msKxmuhKYhiBbUl4JuVJZrbSrDoFpqLRqcYF6XAKoB48PihZQBm8EsgVUdxIgj+FR7mym38OWHT1sXt9Eg1/5MxsuwfYBFfVvIGPRM5EHoDpWxdpaP7GAh+agaskE3C6eGs4QC3hUUgap4oI3BTxEvAYFW7J/LMe3Dp8E2QssRU16BZ1raSRzvI8XVwQAep6dzimqmJ2KoDk7tLYrV9t66OYV8sSeL9IU8GDZkkShxK6Tk6SLKjZG1Oum0e8sJXfyZ6qEyFkYGuBA37hK5pDMk3Fia0d9bgazRSxbsrUvSbpokyioFXgIiFX5ME0DCIPhBWxwecF0E1pWR/sFSzHrg8xMaPmBMnZWjbMaAogG6E7l2TYws6AjkS/RncqrHJoTBXAJbr+wjk2NIUxDYLjUq2BLTJfA6xI0O112W8rp8LO5Upn7e8fp6x4EQhAMMj0WSj1EI8Cjp72X//6RR067XzM/aG1XrrZ1i/6PYPfJSbb89gQbW2vY3TfO6+uDNAVU5h1DKH/bNbUBNalVLAFV4PNAWZLrT9ItBIS97LFs7rPLUEwDfvwtNbSGVbfSsiWkLYaSk1BfRWtHE33dw4CLhrYaNtQFWR8NsG0gxWjehVVVrSaeTqagKggeg84aPxvrgzzSnyJpqUh+Pzk2xr74JLggPpKBSRWX+3Nvv3h6Oboh1MIZQwhMQ60ENFyCpoBn2jc6niviFh7es+Mouf5+VJCoZZA4gWrhXAS4IPGdl7yPH3pzO5zcD4ksovPSOXhymjOhtV2Z2taG/o9gy2Nfw2i6md2941y5KkrIo8YfG/1u9iaytARN9iay7BnNoiaafGD68Ff7aGurmU6evCs+qRI5OB4HTQGTpoCHeK7EkNNqgQk4OUnfZFSJfCJP1OemN1Ng/1iWeL7ExvogDUurVPCnVXXsS2TpGc9NV5hcqcz9fUkeOBSHVB419mmhWic5MKOUpORwKs+ek5PEcyXawmpVYX/GYkO96qrucZIu/FVHHR0RH48Ppp2KMAS0q6Xn5IAaiFXD0Ddf9j6KO/8C7qzhD2KTa+YNre3K1PYZh26EED4hxG4hxD4hxEEhxL84+5cLIXYJIV4QQvyXEMJ09nud7SPO8bbZLcJ8kOIrb2yFVIYNdUEOjOWc+BcW7VU+BrNFvIbAOqkWdBCtgqBJe5WP9dEAN6+opVCWXNca4fKLY+CqByyypTLdqTzr6gJc31ajss7jB/KQHces8kFNgK54hq5DcQZSeax8ke5UnqGsxWDWYn00wOVNYTAE3ak8saBJc8jkxuU1EPCg/H/HUB4DBaDMR65s576j4yQLNq9vCPGBNQ1c0xIhnlO+xt1JFcvkkf4kO0cyZO0yTw6n+dTDh1GBnJwAV8k96jt9bTD03Fney/mrCFrbp0NruxK0/WLOZoy+AFwppVwHrAfeIoS4FPgccJeUchWqRLc7598OjEspVwJ3OedVHHf++G46Lozx+b2Dyh/XWWQxJZx0UUXO27huFXf/6Spu66gDoDPiJxY0yTgLLNqrvLz/TStAhBk6meHJ4Qx96QI3r6jlc5e1wZI6CDUDLqyBIzA+AqMTwCCcTLGpuZpsqczO4TRXLa1iKGsR9hhUm256JvJ0J/OsjvjVmOrQGDPxOVxEl7dy9esupmeiwOVNYd7aUs01LdUADExa7I5nuCwWosHvoWs8x8b6EN+8fDmWLbnv6DjYcZRfdDUQR1WKZta0R4GXSRu3cNDaPg1a2xWh7d/jjIZeKjLOpsd5SeBK4L+d/d8D3u58vt7Zxjl+lRBiXlzbZpdu1TIZzmBLGMwWydpl1tUFCHtcXLcswgfWxWiv8tIzUeCenjEOjGU5MJbFdAnl4mUIkgWb5pDJmtWNRKMBvC6VGzNp2cQCHu5Y3Uj70mqoroXQUvBEgDTK1ew4O383zBsaQ1zdEiHsMYgF1BLulFUiMVmkN1NgJF90EiSEAB/UtEJNCxHTzdXN1by1pZq3t9WQtctOPk7JaL5EyGPQNZ4nUSipa+mspy1s8uGf7mPP/kOohSM1QBQV96OWDWtbOHBw57w9lVeC1vZLobW92LX9Ys7K60YIYQgh9qL+tf0C6AGSUsqSc8oAyvcI570fwDmeQt2tF3/nHUKIZ4UQz766Iswfn3/037jt8uVs2dXPPT0JeiYKWLZKkmAIwYb6IJ0RP//y+FGsgZMw1E26qLLmNIdMetMF+tIFHuhNcvOKWm5pj1Ioy+nMNbaEjfVBbl0V5YaLGiFkQl0Q1TW1AAnlSbafmKAkJf2TFk8Np9k/nlPjlVaJR46nSBZsrm2NsG5ZBJbUQrYIbhcbG4I0+N3YUrmR9UwUeHI4zZ7RScIeddwQ0BTwsLY2gGkItp+YUL9L3rkLNcwsHlmC3+0CnpmPx/FHobV9erS2YbFr+1Re0YIpIUQE+CnwKeA7ThcWIUQL8LCUcq0Q4iDwp1LKAedYD7BRSpl4me9dNItK/pAVwAYgBaFGNq2MckdnA363C8suM5Qt4ne72DaQYs/JLC0hk+vb1BhhrlRWmWxGJ7lhRZTrWiM81JecniwKewwMIUhaJcJOPO5P/3aQ+Oik8kEGNjVXE8+VlL+wZdNc42cgmYPRk6gGqhuE8//c5QK3i+olYVK5Ite31xL1udnalyQ+nIaCDXYeGmq5dkUNTQGTzUtUzk27LDmUzPHFPSdg9DiqQvhRFSEOhLn6dZey7Zl/n+P7f3acaVGJ1vbp0NquBG3DK/Sjl1ImgceBS4GIEGLKa6eZmTxaA0ALgHO8GjVDUqEc5Qt/vhZwQSbOzr3P83/3DeF1sue0hr0kCzZXxKqmExrkSmXaQiaGgLW1AT6yvomkVWJvIsvaWj/dKTX+OLXKr73Kp7LbCKHGGQs25Ev4/R6Slk3PcBqGJiCRYWAoDaPjzIg0BfI4yAGwT0AhQ8oqcU1bDb3pAt/eeZz4sTHI9oE9osphCOyyGmON50oYQoWh3dqXhNGpgE4u1HhoBtUCambbM93z9AxePVrbp0NruxK0DWfndVPvtHYQQviBNwNdwA7gJue024AHnM8POts4xx+TCyHOwizy4Z/+nA1rL0C5m1l0HTrAX//qmPJOsMtsqFerBjfUB7DK0on655teVJKybOq8bnrTqntpCEHE62b/WI6vHhiha1wlXeiZyLM3kYVCATIpcv1DdB06DoUSmG6gDJNJZrwFTOcd1DijCS4/rWEviUKJfcNptZRcZoBGiLRAdYhbVkbZHAspv2dDNRYeH0rTfVSljVNuZhJl44aANi7fsBR4aO5u+jlAa/vMaG0vTm2/mDMO3QghLkZNQBmofwz3Sin/VQixAvgRUAs8B9wqpSwIIXzA94HXoO7Wu6SUR8/wG4u/sgRu5f2blvEf2+9HdftqwNfINRfHuHVVVC3ZNlyO54JLzfYnc8RzJY5M5KfzVF7VVEVnjY/etMUVTWG6xnM8OZzh7W0RWoKmk3YNwh4XpstFPFfElpK/e/o41sALqMdkQGQJpC2wsxCJgNsF+RL43IQCJlctVa2wsMfgyeE0GxuCtARVQKehrEVzyGRoUnXN47ki/3DfHtSQtIGKQ5IBTgIRNq7byO59vwAOztPNPzOn695qbZ8lWtssNm2/GB3U7BzylZv+jb9/6HkoPIeKvNcEwgSXi6/8+WraQl76Jy28LvVcejMFTJcLq6xyZ/ZMFGgJmgw6bmRhj0FnjY/uVN6ZKIINdUHW1voxDRctQZORXJGtx5Pc0xWH0UnU8mwDllSrMK7pApTzzMQSkRD084WrV+F3qyTHU4JX46YqXZrXae08PZLhi9teALkH1Z1tVd/PAGASXf4GEsk8jH93Lm/1K0YHNXt1aG0vXLShP6dMjdu9PDvu+AJvuvs5yO9HdQPbUQL1809/9lo21AUplCUZy8ZwqWBPvekC+xJZABr8HtrCJv0Zi4hXZfyJmAbxXIltAyn6Jy3iuSKxgEnPRB7T5VIrEIs5CIUhMwbBWtYtr2HfASfdmVEFDUFuWF7LfU/34G+J8k+vbcLrUhXBLXAqhHI1s8rqcRwaz/H5R7uBI6ixympmxkf9wCquvbSTrb/5wrm+2eccbehfDq3tSte2NvSzwLN33sUl338OJvaj/IKbnPcG7rrhjcSc+N0AmaLNSK5Ib9oi4jXY1p8i6nPzhsYQbc5kF8C3Dp8kPjiBWaeWbLeFTNLFMkNjWeVS5jVU99VKoMZTDVT2ngmgSv1+dAUkJrjlspVc31bD0KRqXRXKUo2nFkr43S5Mw8X3Do/yyO79qNaNAQScMqiMOtDJrX+yih/86v/M1W19VWhDf27Q2l54aEM/b0T4xs0f5m/ufQw4jGotRFEJkP1cvmEdf3dRAwVbkrRsbKmy0ifyJeK5IntGJ9k/lqMt7KUpoPZv60tCuQzFMu+/ZCkR0026aDOYVcvEDwxOqEh9dgLlJOJHucdZ4AtA3okDsqSW93U2TCeAngr0BJC0bNwCPvSbfjLHu1ELVwSqMoyhWj4GsILWjg76RjKQ+h6LAW3ozxVa2wsNHY9+3kiStGz+6c828+mfeYBe1Ez+UmCYJ/bs5InnV3L39Z2EPS4VFjav0prZUrI+GmB9NIBdhkJZxe9ujfjIFG1CHpVBpy3kxXAJOoo2F1T7aK/y8kDPGIY7hJ1rgmSONasbSRZU6+qqpS3sTWTZN5SejsXdmy5gGmos03SpSH4f29VP5vhvUV1zP6oS96O69m4gxIa1q9Uy+EVSETTnEq3txYhu0c8m1bfxq3ddzJ9840HUrH4zM0u8LWAlH3nLRdOtm+5UnoItaa/yMjBpEfW6ydplIqaBZUsVxxvIlsrYznMzhCBiGvRmCgQM5fFguAQZy2YkX6RrPE/PRJ73r24gVyozWlBZfnbHJ9lQH6Sj2odll/mHh7qgOIwK7RJBtXCmEiCXUS5sMVhysdqMf2tOb+WrRbfozzFa2wsGPXSzIHgnV7+umW3PPIlqRXiBF1CKCgPLWbdmOTcur6W9yktXMkdHtY+kZRP1uhnKFmkOmaQKJdLFMtlSGa8h6KhWC01MQ2BL6E7m6axRC1aiPjfpYplkoUTSsllb6yfkMehLF6aPGQLuPTrG1sOjMD6G6hJ7UWOeXlSr5wTT0fsIgeti3nf5cr6947PzcB9fHdrQzwZa2wsBbegXCM0r/5ZbVkb5/KMPo4RVZiqMqvLbXQp1y9gQq6I5ZPJXHXXTWeetsqQz4iddtKcTI1i2ErxpuAgYLvonLZJWadpXuH/S4g2NISKmQb8zKbWm1s/uuMoYNJCx2D+e48vbnnKuI4Rq1fhRXVgDOIYK7KQCOlH3Rq5ZEeWR3QvfC+F0aEM/O2htzz96jH6BMHDk6zRd/K9El7+OxLGHmVmQXHDeR2BUsmc0xB78WLbkiqYwIY9BxKvydIY9LtJFNaHUXuNnKFskYhrkSiqQVMRU45pWucxlS0KYxkx+z/fsOMp7VkWJmG6+2z3KgYPHUUKvA+9SteBkMo7yOsigWmPjznVGgE7et3bJomztaGYXre3Fgc4ZO0f8w32/5NZVdUAjqiVho6LjTS32SKBaGhke2X2ch/tT5EplZ0KrTMGWTvozNcFlCIg7IV970wUShRLdyTz9GQtbQrJQYs/JSXaOZHh9fZCCLXnoeJIDB4fAG6T9gksgtkTl/ZycQHkeJFEZ7o+iKqoJXMD1mzr49o6FGdBJM/9obS989NDNHNK88m/J2WUSx55BjRG6nFeYmYrhAwRULyFU7ePm9lrW1vjprPFPJ23efmKCrvEcfreLnonCdATAeH9S5dbMTPKBN3dilSXXLosQ8RrTnggP9SW5vzdJKmup1s5QEiX8BCq2RwEVGTAHvI47rlzHlsf2Ag/P/Q07h+ihm9lFa3v+0EM3C4yBI1/nc2//F7aYm+g5/CSqQrhR7mmjqErhLOBIJcikPHz7eBwiVbQ2hHhDY4j10QCWXZ5eRUi6AJlhclQBLij2A0uxypInh9Pc3zvOmlo/fWmLnrEsGII10SC2lGSOn0SJPo+asCqp73Dc5a685AIOjedY7BVBM/tobS9sdIt+ron9NV+4dBkf3tEDye2oFkYRNV4IavIohKokGdQkEsyMsgkwa6HaR6zKx9BACgo5VGslB5EWblsXo2eiwFNDE1AqO+nZioAPPAbYEspplA/xxCm/NZVgqRWz+SJuaa/le08sjtWBZ0K36OcAre15QXvdLFhez93vvp6//NlhSO1CjR9eiBLmBGqs049qBU3FILFQLSKL6eBOeIAA1IRprQ+SLNhOjG8v+0YyhPweMsdHUb7NOWZWAZZQLR1xyqvsfN9FXHtpK00Bky2PfWYO7sXcoA39XKG1PddoQ7+gaeCb7/wg/+/gCPsOHGHGE2AqQ10OaEDF4agCT4hoczVtIS+dNX6SVomtRxIwloPypPpbTwCqfJBIQiAIQRNO9jBTAab0kEO1dFwov+ISsBR/y3Lev7qBe3vGGDjy9bm8GbOONvRzidb2XKIN/YJnE795/zu4p2eML//yCJQPMRODY2pZdgkl2CIqrkgJ1RUOoVopoLqmRWbGId0or4KpbvLzznfmnHMCzt+WgRjUx7htdQMbG0LcuaMHRr8z+0WfY7Shn2u0tucKbegXDdfx09s2053K87GnemH0MHCcmUUnUy0Uj/NeQHVzLZSwTZTIOeWzWx0LroLJY865trN/KTTUcMsFdTQHTcIeg089tAN4Yo7KO/doQz9faG3PNtrQLyqW8qP3/B1el+Anx8b5wVO9UD6IEv5UTI4gqvWSRVWSqcpQQI1RTrV2psY+w+Bpg+IJdX6okSsvrKfO66bB76ElZPL9FxIcOPiVuSzovKAN/XyitT2baEO/KLmWHXdcgV2Gx4cm+PTuATj5AmpsM8uMF8FUS6jsvARqkkuiusMe572a6PJW2qt8bGoMcd2yCHsTWb7/wij7DlR+JZhCG/qFgNb2bKAN/aLnZn72vktIF8sq7+Z4ju5Unm0DKRhOg8yiKogE/OCNQJWX1ho/1y1TLm1tYS9Rn1ou8d5Huxd8WrTZQhv6hYbW9rlCG/qK4ia+cfN6YgETy4njrZaMl9k/nuOyJWHSRZuwx8AuS0KmwSPHk3zviS7ggfm++HlHG/qFjNb2q+GcGnohhAE8C5yQUl4rhFgO/AioBX4LvEdKaQkhvMDdwAbU2uN3Sil7z/DdujK8YmpR/sn1qDHMIipT/dH5vKgFy0tVhtnUtfP9WtuvGK3tV8LZGPpXEtTsg0DXKdufA+6SUq5COcre7uy/HRiXUq4E7nLO05xzxoCnUS2a+4CH0BXhj0LresGhtX3OkVKe8YVKH7MduBLYipodGQXczvFNwM+dzz8HNjmf3c554gzfL/VLv2bzNR+61trWr7l4nY0NP9sW/ZeAjzLj2xQFklLKKQfXAVTSSJz3ftQVTC2Fi57l72g0c4nWtea84IyGXghxLRCXUu45dfdpTpVncezU771DCPGsEOLZs7pSjeYcMlu6dr5ba1uzoDibMMVvBN4mhHgraiVDFaolFBFCuJ3WTTMqFiioVlALMCCEcKPWNo+9+EullFuABzY5LAAABZpJREFULaAnrDTzwqzoGrS2NQuPM7bopZSfkFI2SynbgHcBj0kp/wLYAdzknHYbM35ODzrbOMcfkwvBh1OjOQWta835xKtJJfgx4B+FEEdQY5XfcvZ/C4g6+/8R+Piru0SNZk7RutZUHHrBlOa8QC+Y0lQq59qPXqPRaDSLEG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqHG3oNRqNpsLRhl6j0WgqnLMy9EKIXiHEfiHEXiHEs86+WiHEL4QQLzjvNc5+IYT4ihDiiBDid0KI185mATSaV4PWtuZ84JW06N8kpVwvpbzE2f44sF1KuQrY7mwDXAOscl53AP9xri5Wo5kltLY1Fc2rGbq5Hvie8/l7wNtP2X+3VPwGiAghYq/idzSauUZrW1NRnK2hl8A2IcQeIcQdzr5GKeUQgPPe4OxfCvSf8rcDzr7fQwhxhxDi2anuskYzT2htayoe91me90Yp5aAQogH4hRDi+Zc5V5xmn/yDHVJuAbYACCH+4LhGM0dobWsqnrNq0UspB533OPBTYCMwMtVtdd7jzukDQMspf94MDJ6rC9ZoziVa25rzgTMaeiFEUAgRnvoMXA0cAB4EbnNOuw14wPn8IPCXjofCpUBqqhus0SwktLY15wtnM3TTCPxUCDF1/g+llI8KIZ4B7hVC3A4cB97hnP8w8FbgCJAF3nsWv5EBDr/Ca1/M1AGj830Rc8RCKGvrS+zX2j63LIRnPZcshPK+lLZ/DyHl/A8hCiGePcW1reI5n8p7PpX1dJxP5T+fygqLq7x6ZaxGo9FUONrQazQaTYWzUAz9lvm+gDnmfCrv+VTW03E+lf98KissovIuiDF6jUaj0cweC6VFr9FoNJpZYt4NvRDiLUKIw05EwI+f+S8WNkKIFiHEDiFElxDioBDig87+io2IKIQwhBDPCSG2OtvLhRC7nLL+lxDCdPZ7ne0jzvG2+bzu2aTSdA1a2872otT2vBp6IYQBfA0VFXA1cIsQYvV8XtM5oAR8SErZCVwK3OmUqZIjIn4Q6Dpl+3PAXU5Zx4Hbnf23A+NSypXAXc55FUeF6hq0tmGxaltKOW8vYBPw81O2PwF8Yj6vaRbK+ADwP1CLZmLOvhhw2Pn8DeCWU86fPm8xvFBhALYDVwJbUfFgRgH3i58x8HNgk/PZ7Zwn5rsMs3BPKl7XTrm0theJtud76OasogEuVpzu22uAXbzKiIgLmC8BHwXKznYUSEopS872qeWZLqtzPOWcX2ks9md6RrS2gUWk7fk29GcVDXAxIoQIAT8B/peUcuLlTj3NvkVxD4QQ1wJxKeWeU3ef5lR5FscqiYoup9b277EotH22YYpni4qMBiiE8KAqwn9KKe9zdo8IIWJSyqEKioj4RuBtQoi3Aj6gCtUKiggh3E7L5tTyTJV1QAjhBqqBsbm/7FlnMT/Tl0Vre3Fqe75b9M8Aq5yZbBN4FypC4KJFqAhZ3wK6pJRfPOVQxUVElFJ+QkrZLKVsQz27x6SUfwHsAG5yTntxWafuwU3O+Qum1XMOqThdg9b2otb2fE8SoKIBdgM9wP+e7+s5B+XZjOqy/Q7Y67zeihqv2w684LzXOucLlIdGD7AfuGS+y/BHlvsKYKvzeQWwGxXl8ceA19nvc7aPOMdXzPd1z+L9qChdO2XS2l6k2tYrYzUajabCme+hG41Go9HMMtrQazQaTYWjDb1Go9FUONrQazQaTYWjDb1Go9FUONrQazQaTYWjDb1Go9FUONrQazQaTYXz/wFi6orcLiFyXAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# example of a \"bad data point\" (i.e. (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0) and (int(dcm.RescaleIntercept) > -100) == True)\n",
"example_img = train_img_path + train_df.index[102] + \".dcm\"\n",
"\n",
"dicom = pydicom.dcmread(example_img)\n",
"\n",
"fig, ax = plt.subplots(1, 2)\n",
"\n",
"ax[0].imshow(window_testing(dicom, window_without_correction), cmap=plt.cm.bone);\n",
"ax[0].set_title(\"original\")\n",
"ax[1].imshow(window_testing(dicom, window_with_correction), cmap=plt.cm.bone);\n",
"ax[1].set_title(\"corrected\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load Data"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"_cell_guid": "79c7e3d0-c299-4dcb-8224-4455121ee9b0",
"_uuid": "d629ff2d2480ee46fbb7e2d37f6b5fab8052498a"
},
"outputs": [],
"source": [
"def _read(path, desired_size):\n",
" \"\"\"Will be used in DataGenerator\"\"\"\n",
" \n",
" dcm = pydicom.dcmread(path)\n",
" \n",
" try:\n",
" img = bsb_window(dcm)\n",
" except:\n",
" img = np.zeros(desired_size)\n",
" \n",
" \n",
" img = cv2.resize(img, desired_size[:2], interpolation=cv2.INTER_LINEAR)\n",
" \n",
" return img\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Data generators\n",
"\n",
"Inherits from keras.utils.Sequence object and thus should be safe for multiprocessing.\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"class DataGenerator(keras.utils.Sequence):\n",
"\n",
" def __init__(self, list_IDs, labels=None, batch_size=1, img_size=(512, 512, 1), \n",
" img_dir=train_img_path, *args, **kwargs):\n",
"\n",
" self.list_IDs = list_IDs\n",
" self.labels = labels\n",
" self.batch_size = batch_size\n",
" self.img_size = img_size\n",
" self.img_dir = img_dir\n",
" self.on_epoch_end()\n",
"\n",
" def __len__(self):\n",
" return int(ceil(len(self.indices) / self.batch_size))\n",
"\n",
" def __getitem__(self, index):\n",
" indices = self.indices[index*self.batch_size:(index+1)*self.batch_size]\n",
" list_IDs_temp = [self.list_IDs[k] for k in indices]\n",
" \n",
" if self.labels is not None:\n",
" X, Y = self.__data_generation(list_IDs_temp)\n",
" return X, Y\n",
" else:\n",
" X = self.__data_generation(list_IDs_temp)\n",
" return X\n",
" \n",
" def on_epoch_end(self):\n",
" \n",
" \n",
" if self.labels is not None: # for training phase we undersample and shuffle\n",
" # keep probability of any=0 and any=1\n",
" keep_prob = self.labels.iloc[:, 0].map({0: 0.35, 1: 0.5})\n",
" keep = (keep_prob > np.random.rand(len(keep_prob)))\n",
" self.indices = np.arange(len(self.list_IDs))[keep]\n",
" np.random.shuffle(self.indices)\n",
" else:\n",
" self.indices = np.arange(len(self.list_IDs))\n",
"\n",
" def __data_generation(self, list_IDs_temp):\n",
" X = np.empty((self.batch_size, *self.img_size))\n",
" \n",
" if self.labels is not None: # training phase\n",
" Y = np.empty((self.batch_size, 6), dtype=np.float32)\n",
" \n",
" for i, ID in enumerate(list_IDs_temp):\n",
" X[i,] = _read(self.img_dir+ID+\".dcm\", self.img_size)\n",
" Y[i,] = self.labels.loc[ID].values\n",
" \n",
" return X, Y\n",
" \n",
" else: # test phase\n",
" for i, ID in enumerate(list_IDs_temp):\n",
" X[i,] = _read(self.img_dir+ID+\".dcm\", self.img_size)\n",
" \n",
" return X"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Metrics"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"from keras import backend as K\n",
"\n",
"def weighted_log_loss(y_true, y_pred):\n",
" \"\"\"\n",
" Can be used as the loss function in model.compile()\n",
" ---------------------------------------------------\n",
" \"\"\"\n",
" \n",
" class_weights = np.array([2., 1., 1., 1., 1., 1.])\n",
" \n",
" eps = K.epsilon()\n",
" \n",
" y_pred = K.clip(y_pred, eps, 1.0-eps)\n",
"\n",
" out = -( y_true * K.log( y_pred) * class_weights\n",
" + (1.0 - y_true) * K.log(1.0 - y_pred) * class_weights)\n",
" \n",
" return K.mean(out, axis=-1)\n",
"\n",
"\n",
"def _normalized_weighted_average(arr, weights=None):\n",
" \"\"\"\n",
" A simple Keras implementation that mimics that of \n",
" numpy.average(), specifically for this competition\n",
" \"\"\"\n",
" \n",
" if weights is not None:\n",
" scl = K.sum(weights)\n",
" weights = K.expand_dims(weights, axis=1)\n",
" return K.sum(K.dot(arr, weights), axis=1) / scl\n",
" return K.mean(arr, axis=1)\n",
"\n",
"\n",
"def weighted_loss(y_true, y_pred):\n",
" \"\"\"\n",
" Will be used as the metric in model.compile()\n",
" ---------------------------------------------\n",
" \n",
" Similar to the custom loss function 'weighted_log_loss()' above\n",
" but with normalized weights, which should be very similar \n",
" to the official competition metric:\n",
" https://www.kaggle.com/kambarakun/lb-probe-weights-n-of-positives-scoring\n",
" and hence:\n",
" sklearn.metrics.log_loss with sample weights\n",
" \"\"\"\n",
" \n",
" class_weights = K.variable([2., 1., 1., 1., 1., 1.])\n",
" \n",
" eps = K.epsilon()\n",
" \n",
" y_pred = K.clip(y_pred, eps, 1.0-eps)\n",
"\n",
" loss = -( y_true * K.log( y_pred)\n",
" + (1.0 - y_true) * K.log(1.0 - y_pred))\n",
" \n",
" loss_samples = _normalized_weighted_average(loss, class_weights)\n",
" \n",
" return K.mean(loss_samples)\n",
"\n",
"\n",
"def weighted_log_loss_metric(trues, preds):\n",
" \"\"\"\n",
" Will be used to calculate the log loss \n",
" of the validation set in PredictionCheckpoint()\n",
" ------------------------------------------\n",
" \"\"\"\n",
" class_weights = [2., 1., 1., 1., 1., 1.]\n",
" \n",
" epsilon = 1e-7\n",
" \n",
" preds = np.clip(preds, epsilon, 1-epsilon)\n",
" loss = trues * np.log(preds) + (1 - trues) * np.log(1 - preds)\n",
" loss_samples = np.average(loss, axis=1, weights=class_weights)\n",
"\n",
" return - loss_samples.mean()\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"\n",
"class PredictionCheckpoint(keras.callbacks.Callback):\n",
" \n",
" def __init__(self, test_df, valid_df, \n",
" test_images_dir=test_img_path, \n",
" valid_images_dir=train_img_path, \n",
" batch_size=32, input_size=(224, 224, 3)):\n",
" \n",
" self.test_df = test_df\n",
" self.valid_df = valid_df\n",
" self.test_images_dir = test_images_dir\n",
" self.valid_images_dir = valid_images_dir\n",
" self.batch_size = batch_size\n",
" self.input_size = input_size\n",
" \n",
" def on_train_begin(self, logs={}):\n",
" self.test_predictions = []\n",
" self.valid_predictions = []\n",
" \n",
" def on_epoch_end(self,batch, logs={}):\n",
" self.test_predictions.append(\n",
" self.model.predict_generator(\n",
" DataGenerator(self.test_df.index, None, self.batch_size, self.input_size, self.test_images_dir), verbose=2)[:len(self.test_df)])\n",
" \n",
"\n",
"\n",
"class ResNet50_Model:\n",
" \n",
" def __init__(self, engine, input_dims, batch_size=5, num_epochs=4, learning_rate=1e-3, \n",
" decay_rate=1e-6, decay_steps=1, weights=\"imagenet\", verbose=1):\n",
" \n",
" self.engine = engine\n",
" self.input_dims = input_dims\n",
" self.batch_size = batch_size\n",
" self.num_epochs = num_epochs\n",
" self.learning_rate = learning_rate\n",
" self.decay_rate = decay_rate\n",
" self.decay_steps = decay_steps\n",
" self.weights = weights\n",
" self.verbose = verbose\n",
" self._build()\n",
"\n",
" def _build(self):\n",
" \n",
" \n",
" engine = self.engine(include_top=False, weights=self.weights, input_shape=self.input_dims,\n",
" backend = keras.backend, layers = keras.layers,\n",
" models = keras.models, utils = keras.utils)\n",
" \n",
" x = keras.layers.GlobalAveragePooling2D(name='avg_pool')(engine.output)\n",
" out = keras.layers.Dense(6, activation=\"sigmoid\", name='dense_output')(x)\n",
"\n",
" self.model = keras.models.Model(inputs=engine.input, outputs=out)\n",
"\n",
" self.model.compile(loss=\"binary_crossentropy\", optimizer=keras.optimizers.Adam(), metrics=[weighted_loss])\n",
" \n",
"\n",
" def fit_and_predict(self, train_df, valid_df, test_df):\n",
" \n",
" # callbacks\n",
" pred_history = PredictionCheckpoint(test_df, valid_df, input_size=self.input_dims)\n",
" scheduler = keras.callbacks.LearningRateScheduler(lambda epoch: self.learning_rate * pow(self.decay_rate, floor(epoch / self.decay_steps)))\n",
" \n",
" self.model.fit_generator(\n",
" DataGenerator(\n",
" train_df.index, \n",
" train_df, \n",
" self.batch_size, \n",
" self.input_dims, \n",
" train_img_path\n",
" ),\n",
" epochs=self.num_epochs,\n",
" verbose=self.verbose,\n",
" use_multiprocessing=True,\n",
" workers=4,\n",
" callbacks=[pred_history, scheduler]\n",
" )\n",
" \n",
" return pred_history\n",
" \n",
" def save(self, path):\n",
" self.model.save_weights(path)\n",
" \n",
" def load(self, path):\n",
" self.model.load_weights(path)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Train model and predict\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Epoch 1/5\n",
" 61/7872 [..............................] - ETA: 1:52:11 - loss: 0.2221 - weighted_loss: 0.2492"
]
}
],
"source": [
"# Train/Valid/Test Split\n",
"train_df_ss = ShuffleSplit(n_splits=10, test_size=0.1, random_state=42).split(train_df.index)\n",
"\n",
"\n",
"train_idx, valid_idx = next(train_df_ss)\n",
"\n",
"# model\n",
"model = ResNet50_Model(engine=ResNet50, input_dims=(224, 224, 3), batch_size=32, learning_rate=5e-4,\n",
" num_epochs=5, decay_rate=0.8, decay_steps=1, weights=\"imagenet\", verbose=1)\n",
"\n",
"#predictions \n",
"history = model.fit_and_predict(train_df.iloc[train_idx], train_df.iloc[valid_idx], test_df)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Ensemble and average all submission_predictions."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"test_df.iloc[:, :] = np.average(history.test_predictions, axis=0, weights=[0, 1, 2, 4, 6]) # let's do a weighted average for epochs (>1)\n",
"\n",
"test_df = test_df.stack().reset_index()\n",
"\n",
"test_df.insert(loc=0, column='ID', value=test_df['Image'].astype(str) + \"_\" + test_df['Diagnosis'])\n",
"\n",
"test_df = test_df.drop([\"Image\", \"Diagnosis\"], axis=1)\n",
"\n",
"test_df.to_csv('submission.csv', index=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from IPython.display import FileLink, FileLinks\n",
"FileLink('submission.csv')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}